cleanup before you start

rm( list=ls() )
library(moments)

setup

set.seed(1234)
N  <- 500
H  <- 5 

phi0    <- 0.2
phi1    <- 0.97
theta1  <- 0.6
sig.eps <- 0.1

simulate DGP

y  <- rep(0,N+H)
eps     <- rnorm(N+H,0,sig.eps)

y[1] <- phi0/(1.0-phi1) + eps[1]

for( t in 2:(N+H) ){
    y[t] <- phi0 + phi1*y[t-1] + eps[t] - theta1*eps[t-1]
}

y.out <- y[(N+1):(N+H)]
y     <- y[1:N]

plot( y , t='p' , pch=16, col='darkorange2' , tck = 0.02 , xlim=c(1,N+H) )
grid( lwd=1 , col="darkgrey" )

LBQ Test

Box.test( y, lag=22 , type="Ljung-Box" )

## 
##  Box-Ljung test
## 
## data:  y
## X-squared = 4597.3, df = 22, p-value < 2.2e-16

# ACF & PACF
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( y , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( y , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)

Estimate ARMA Models

ar1    = arima(y,order=c(1,0,0))
ar2    = arima(y,order=c(2,0,0))
ma1    = arima(y,order=c(0,0,1))
arma11 = arima(y,order=c(1,0,1))

ar1_aic    <- (-2*ar1$loglik+2*3)/N
ar2_aic    <- (-2*ar2$loglik+2*4)/N
ma1_aic    <- (-2*ma1$loglik+2*3)/N
arma11_aic <- (-2*arma11$loglik+2*4)/N

ar1_bic    <- (-2*ar1$loglik+log(N)*3)/N
ar2_bic    <- (-2*ar2$loglik+log(N)*4)/N
ma1_bic    <- (-2*ma1$loglik+log(N)*3)/N
arma11_bic <- (-2*arma11$loglik+log(N)*4)/N

round( rbind( c(ar1$loglik,ar2$loglik,ma1$loglik,arma11$loglik), 
              c(ar1_aic,ar2_aic,ma1_aic,arma11_aic) , 
              c(ar1_bic,ar2_bic,ma1_bic,arma11_bic) ) ,  3 )

##         [,1]    [,2]    [,3]    [,4]
## [1,] 383.158 410.141 234.417 425.704
## [2,]  -1.521  -1.625  -0.926  -1.687
## [3,]  -1.495  -1.591  -0.900  -1.653

FITTED VALUES

ar1_mu     <- y-ar1$residuals
ar2_mu     <- y-ar2$residuals
ma1_mu     <- y-ma1$residuals
arma11_mu  <- y-arma11$residuals
ar1_res    <- as.numeric(ar1$residuals)
ar2_res    <- as.numeric(ar2$residuals)
ma1_res    <- as.numeric(ma1$residuals)
arma11_res <- as.numeric(arma11$residuals)

Plot a bunch of stuff..

# AR1
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
plot( y , t='p' , pch=16, col='darkorange2' , tck = 0.02 , xlim=c(1,N+H) )
lines( ar1_mu , t='l' , lwd=2 , col='blue3' )
grid( lwd=1 , col="darkgrey" )
plot( ar1_res , col='purple' )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )

par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( ar1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( ar1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)

# AR2
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
plot( y , t='p' , pch=16, col='darkorange2' , tck = 0.02 , xlim=c(1,N+H) )
lines( ar2_mu , t='l' , lwd=2 , col='blue3' )
grid( lwd=1 , col="darkgrey" )
plot( ar2_res , col='purple' )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )

par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( ar2_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( ar2_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)

# MA1
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
plot( y , t='p' , pch=16, col='darkorange2' , tck = 0.02 , xlim=c(1,N+H) )
lines( ma1_mu , t='l' , lwd=2 , col='blue3' )
grid( lwd=1 , col="darkgrey" )
plot( ma1_res , col='purple' )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )

par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( ma1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( ma1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)

# ARMA11
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
plot( y , t='p' , pch=16, col='darkorange2' , tck = 0.02 , xlim=c(1,N+H) )
lines( arma11_mu , t='l' , lwd=2 , col='blue3' )
grid( lwd=1 , col="darkgrey" )
plot( arma11_res , col='purple' )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )

par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( arma11_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( arma11_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)

Residuals

par( mar=c(2,2,3,2) , mfrow=c(2,2) )
kernel <- density(ar1_res/sqrt(ar1$sigma2))
plot( kernel , main='AR1' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )

kernel <- density(ar2_res/sqrt(ar2$sigma2))
plot( kernel , main='AR2' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )

kernel <- density(ma1_res/sqrt(ma1$sigma2))
plot( kernel , main='MA1' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )

kernel <- density(arma11_res/sqrt(arma11$sigma2))
plot( kernel , main='ARMA11' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )

par( mar=c(2,2,3,2) , mfrow=c(2,2) )

qqnorm(ar1_res,col='tomato',main='AR1')
qqline(ar1_res,lwd=2,lty=3)
qqnorm(ar2_res,col='tomato',main='AR2')
qqline(ar2_res,lwd=2,lty=3)
qqnorm(ma1_res,col='tomato',main='MA1')
qqline(ma1_res,lwd=2,lty=3)
qqnorm(arma11_res,col='tomato',main='ARMA11')
qqline(arma11_res,lwd=2,lty=3)

Residual Test

ar1.res.jb    <- jarque.test(ar1_res)$p.value
ar2.res.jb    <- jarque.test(ar2_res)$p.value
ma1.res.jb    <- jarque.test(ma1_res)$p.value
arma11.res.jb <- jarque.test(arma11_res)$p.value

ar1.res.lbq    <- Box.test( ar1_res, lag=22 , type="Ljung-Box" )$p.value
ar2.res.lbq    <- Box.test( ar2_res, lag=22 , type="Ljung-Box" )$p.value
ma1.res.lbq    <- Box.test( ma1_res, lag=22 , type="Ljung-Box" )$p.value
arma11.res.lbq <- Box.test( arma11_res, lag=22 , type="Ljung-Box" )$p.value

round( rbind( c(ar1.res.jb,ar2.res.jb,ma1.res.jb,arma11.res.jb) , 
       c(ar1.res.lbq,ar2.res.lbq,ma1.res.lbq,arma11.res.lbq) ) , 3 )

##       [,1]  [,2]  [,3]  [,4]
## [1,] 0.123 0.436 0.001 0.424
## [2,] 0.000 0.000 0.000 0.452

FORECAST

ar1_pred    <- predict( ar1 , n.ahead=H )
ar2_pred    <- predict( ar2 , n.ahead=H )
ma1_pred    <- predict( ma1 , n.ahead=H )
arma11_pred <- predict( arma11 , n.ahead=H )

ar1_mse    = mean( (y.out - as.numeric(ar1_pred$pred) )**2 )
ar2_mse    = mean( (y.out - as.numeric(ar2_pred$pred) )**2 )
ma1_mse    = mean( (y.out - as.numeric(ma1_pred$pred) )**2 )
arma11_mse = mean( (y.out - as.numeric(arma11_pred$pred) )**2 )

plot( c((N-10):(N+H)) , c(y[(N-10):N], y.out) , main=sprintf('AR(1) MSE %3.3f',ar1_mse) , ylim=c(min(y),max(y)) , ylab='',xlab='', tck = 0.02 , pch=16 , col='darkorange' ) 
abline( v=N , lwd=2 )
abline( h=ar1$coef['intercept'] , lwd=2 )
grid( lwd=1 , col="darkgrey" )
lines( c((N-10):N) , ar1_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( c((N+1):(N+H)) , as.numeric(ar1_pred$pred)  , t='b' , lwd=2 , col='blue3' )

plot( c((N-10):(N+H)) , c(y[(N-10):N], y.out) , main=sprintf('AR(2) MSE %3.3f',ar2_mse) , ylim=c(min(y),max(y)) , ylab='',xlab='', tck = 0.02 , pch=16 , col='darkorange') 
abline( v=N , lwd=2 )
abline( h=ar2$coef['intercept'] , lwd=2 )
grid( lwd=1 , col="darkgrey" )
lines( c((N-10):N) , ar2_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( c((N+1):(N+H)) , as.numeric(ar2_pred$pred)  , t='b' , lwd=2 , col='blue3' )

plot( c((N-10):(N+H)) , c(y[(N-10):N], y.out) , main=sprintf('MA(1) MSE %3.3f',ma1_mse) , ylim=c(min(y),max(y)) , ylab='',xlab='', tck = 0.02 , pch=16 , col='darkorange') 
abline( v=N , lwd=2 )
abline( h=ma1$coef['intercept'] , lwd=2 )
grid( lwd=1 , col="darkgrey" )
lines( c((N-10):N) , ma1_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( c((N+1):(N+H)) , as.numeric(ma1_pred$pred)  , t='b' , lwd=2 , col='blue3' )

plot( c((N-10):(N+H)) , c(y[(N-10):N], y.out) , main=sprintf('ARMA(1,1) MSE %3.3f',arma11_mse) , ylim=c(min(y),max(y)) , ylab='',xlab='', tck = 0.02 , pch=16  , col='darkorange') 
abline( v=N , lwd=2 )
abline( h=ma1$coef['intercept'] , lwd=2 )
grid( lwd=1 , col="darkgrey" )
lines( c((N-10):N) , arma11_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( c((N+1):(N+H)) , as.numeric(arma11_pred$pred)  , t='b' , lwd=2 , col='blue3' )

round(c( ar1_mse , ar2_mse , ma1_mse ,  arma11_mse )*100 , 3 )

## [1] 1.381 1.104 4.084 1.043